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An introduction to linear stability analysis for deciphering spa¬ 
tial patterns in signaling networks 

Jasmine A. Nirody “ and Padmini Rangamani*^ 


Mathematical modeling is now used commonly in the analysis of signaling networks. With advances in 
high resolution microscopy, the spatial location of different signaling molecules and the spatio-temporal 
dynamics of signaling microdomains are now widely acknowledged as key features of biochemical signal 
transduction. Reaction-diffusion mechanisms are commonly used to model such features, often with a heavy 
reliance on numerical simulations to obtain results. However, simulations are parameter dependent and may 
not be able to provide an understanding of the full range of the system responses. Analytical approaches on 
the other hand provide a framework to study the entire phase space. In this tutorial, we provide a largely 
analytical method for studying reaction-diffusion models and analyzing their stability properties. Using two 
representative biological examples, we demonstrate how this approach can guide experimental design. In 
order to make such analyses more accessible, we provide an easy-to-use graphical user interface to test the 
stability behavior of biological systems (available at http://www.ocf.berkley.edu/~jnirody 


1 


1 Introduction 


Q 

Mathematical models have proven to be extremely useful in elucidating the properties of biological sig¬ 
naling networks IS [1210 . A common method of modeling such systems is to use ordinary differential 
equations (ODEs) to describe the dynamics of various signaling components. However, given the compart- 
mental nature of cells, we need to integrate both the spatial and temporal dynamics of signaling in cells. 
Partial differential equation (PDE) models allow us to consider the spatial aspects of biological signaling 
networks. Reaction-diffusion models, in particular, are ubiquitous in biology, particularly in systems which 
give rise to spatial patterning ifTOlfTTlfTSlfTSll . 

PDE models often exhibit a rich behavior as parameters are varied, making analysis of these systems 
challenging. Despite their seeming simplicity, PDE models prove difficult to analyze due to their rich be¬ 
havior across parameter space, especially in comparison to the phase-space of their reaction-only ODE 
counterparts. This daunting complexity often leads to these models being studied primarily through sim¬ 
ulations, which provide instant but transient gratification: while we gain insight into the behavior of these 
systems, the variety of these behaviors across parameter space nearly assures that we cannot grasp their full 
capabilities through simulation alone. 

Several computational tools have arisen over the past years to simulate complex biological models ||T] 
|4l|6l|23l. In the category of analytical methods, nonlinear analysis provides a comprehensive global picture 
of the system, but can be very difficult to implement even for simple models. A recently proposed method, 
local perturbation analysis (EPA), gives information about stimulus-driven cellular responses (we point the 
interested reader to [8] for more details on this method). However, linear (Turing) stability analysis remains 
the gold standard for understanding spontaneous pattern formation in reaction-diffusion systems. Despite 
its well-accepted usefulness in exploring such systems, there are few user-friendly tools available to the 
larger biological community to utilize these analytical methods. In the following, we provide an analytic 
protocol for the study of reaction-diffusion systems. We demonstrate this protocol using two simple, classic 
biological examples of pattern formation and cell polarization. We also note that the method outlined is 
easily generalizable to any signaling network. 

2 General model framework 

2.1 Reaction-diffusion equations 

We consider chemical reactions occurring along a one-dimensional segment (0 < x < L), which serves 
as an approximation of a transection of the cell. We focus primarily on examples corresponding to two 
chemical reactants (Eigure[^), but note that the method described (as well as our provided computational 
tool) is easily extensible to systems with an arbitrary number of species. Depending on the relative values 
of the diffusion coefficient, this extension can serve to model several cytoplasmic species (Eigure[Tp) or a 
mixture of cytoplasmic and membrane-bound species (Eigure[Tp). Here, we provide a general example with 
two species U and A. We denote the concentrations of these species as u and a respectively, both in units of 
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Figure 1: Schematic of a simplified signaling model. (A) Purely nonspatial model corresponds to a well- 
mixed system. The model consists of an interconversion between inactive (U) and active (A) forms of a 
single protein. The active form effects a positive feedback on its own production; many signaling networks 
take the form of a cooperative Hill function. An extended one-dimensional spatial model can be represen¬ 
tative of (B) two cytoplasmic species or (C) one cytoplasmic and one membrane-bound species, depending 
on the relative diffusion coefficients of the forms. 


molecules/length. The partial differential equations are: 


da d'^a 
du ^ d‘^u , . 


( 1 ) 


The choice of functions /(a, u) and g{a, u) is arbitrary. We will first outline the general analysis method 
and then apply it to two specific systems to demonstrate its usefulness. 


2.2 Stability analysis for the well-mixed system 

The first step is to find steady state solutions {a*,u*) of the purely kinetic system (in the case of homogenous 
spatial conditions), we solve for (/(a, u),g{a, u)) = (0, 0). Additionally, to enforce mass conservation, we 
impose the constraint p = uo + uq = u{t) + a{t), where p is a constant denoting the total amount of protein. 
Depending on the form of /(o, u) and g{a, u), multiple steady states can exist. 

These homogenous equilibria also serve as steady states for the spatially extended system. We are 
interested in their stability throughout the parameter space. Next, for the spatially homogenous system, we 
look at the eigenvalues of the Jacobian matrix to analyze the linear stability of these states under small 
perturbations from equilibrium ||26]| . The Jacobian matrix for the purely kinetic system is 

j ^ ( f{a,u)a f{a,u)u \ ^ 2 ) 

V 9{a,u)a g{a,u)u ) ' 

Here, fa and /„ denote the partial derivatives of /(a, u) with respect to a and u, respectively. The lin¬ 
earization of the reaction terms around an equilibrium point {a*,u*) is simply J\(^a*,u*)[da du]'^, where 
J\{a*,u*) represents the Jacobian matrix evaluated at (a*, u*). An equilibrium point is said to be stable if the 
eigenvalues of the J\{a*,u*) have real parts less than or equal to zero and unstable if any one eigenvalue 
has a real part greater than zero. For example, in bistable systems such as MAP kinase, there are two stable 
steady states and one unstable steady state. 
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2.3 Stability analysis for the spatial model 


What happens when we consider a spatially heterogeneous system? In order to obtain the spatially heteroge¬ 
nous solutions, we linearize Equations ([T]l, we arrive at 


d 'da 
dt du 


Da 

0 


0 

Du 


9^ da 
dx^ du 


+ J 


da 

du 


( 3 ) 


Here, J is the Jacobian for the (spatially homogenous) reaction equations. In order to analyze the stability of 
this system with respect to perturbations, we must first note that now these perturbations depend both on time 
and space. A convenient form for such perturbations is [da du] = [da* du*]e~^^e^^^, where the term 
is a common way of representing a spatial wave, with k as the wavenumber (we refer the interested reader 
to Ehll l. Substituting this into the linearized Equation Q leads to the Jacobian of the spatially extended 
system 

f{a,u)a-DAk'^ f{a,u)u 
g{a,u)a g{a,u)u- Duk'^ 

The eigenvalues of this matrix allow us to study the stability properties of the system: nonnegative eigen¬ 
values correspond to a loss of stability of the system. Since the eigenvalue expressions contain the unknown 
variable k, we must consider a one parameter family of solutions, one for each wavenumber. Note that the 
eigenvalues for the well-mixed system are achieved when k = 0. Primarily, we are interested in the emer¬ 
gence of heterogenous patterns that occur via perturbations within a finite range of critical wavenumbers 

0 ^ fee ^ femax- 

The value of these critical wavenumbers becomes important when we are faced with a finite domain 
size. Since cells are of finite size, we want to investigate how the length of the domain affects the system 
response. This is because the spatial pattern that emerges from an instability corresponding to a wavenumber 
fe has wavelength w = ^; accordingly, for finite systems, only values of fe above a certain threshold will 
generate any meaningful spatial patterns. 

We emphasize that the steps discussed in this article are highly generalizable, and we aim to highlight 
how such comprehensive treatments of mathematical models in biology are powerful and insightful. 




3 Examples 

In this section, we demonstrate the above analytical protocol using two examples, each representing a dis¬ 
tinct origin of spatial patterning. Eirst, we discuss Turing’s “diffusion-driven instability” using morphogen¬ 
esis as a classical example. Second, we consider a situation in which spatial patterns arise even though 
the criteria for Turing instabilities are not met, via the “pinning” of a traveling wave. To illustrate this 
phenomena, we turn to a simple model of eukaryotic cell polarization. 

3.1 Classic Turing spots 

Morphogenesis, the process by which an embryo gains its structure during development, is one of the most 
studied problems in biology. In particular, Turing’s model of pattern formation, originally applied to the 
development of animal coat patterns ||29ll . proved useful for describing this phenomenon. Murray llT9ll used 
and analyzed the following system of non-linear reactions, which has since become the classical example of 
the Turing instability. In this system, f{a, u) and g{a, u) are described as follows: 

f{a,u) =b-a- (5) 

g{a,u) =aic-u)- 

where b, p, K, and c are reaction rate parameters. 
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3.1.1 Well-mixed model 

We first begin with the (relatively) easier task of analyzing the well-mixed model before considering the full 
reaction-diffusion system. This system has two homogenous steady states, which are obtained by solving 
/(a, u) = 0 and g{a, u) = 0. We denote them as (a^^, and (o^, We do not explicitly write out the 
expressions for the sake of brevity. 

Next, we construct the Jacobian of this system to obtain 

Ka^pu pu 1 —ap 

(l+a+Ka^)'^ ~ 1+a+Ka^ ~ ^ l+a+Ka'^ 

Ko?pu pu ap 

{l+a+Ka^)'^ l+a+Ka^ “ l+a+Ka^ 

The eigenvalues of this matrix have negative real parts indicating that the purely kinetic system is stable. 

3.1.2 Spatial patterns 

The criteria for an equilibrium point to undergo a diffusion-driven instability are as follows. 

• First, the equilibrium point must be stable for the purely kinetic system and the loss of stability in the 
steady state must be purely spatially dependent. 

• Second, the Jacobian of the spatially extended system must either have a trace less than 0 or a deter¬ 
minant greater than 0. That is, the system cannot undergo a classical Turing instability if both: 

= /a + 5'«<0, (7) 

det(J*) = faQu - fuQa > 0. 

For the above system, the conditions for the Turing mechanism are met, and spatial patterning through 
a diffusion-driven instability can occur. Note that there are six parameters in the system, including 
D = Da/Du, b, c, K, a, and p. In Figure]^ we fix K, a, and p and study the effect of D on the phase 
space admitting Turing patterns over the ranges of parameters b and c. When D = 10, the range of b 
and c over which Turing patterns are observed is small. However, as D increases, the range of b and 
c increases. Biologically, this means when A diffuses much faster than U, the range of parameters b 
and c over which Turing patterns can be observed increases. For more details on linear instabilities 
via the Turing mechanism, we refer the interested reader to lfT9l . 

3.1.3 Key takeaway points 

Classical Turing systems consist of two components, usually an activator and an inhibitor . Since Turing’s 
original paper on pattern formation via a reaction-diffusion mechanism, several models have been proposed 
that are applicable to a wide variety of biochemical problems . 

Above, we have described and analyzed Thomas’ model for substrate inhibition based on experiments 
in which an immobilized enzyme on a membrane react with diffusing substrate and co-substrate molecules 
IITTI . This model, and several others proposed after it El |24l help us understand biologically plausible 
mechanisms of signaling. Some of the key points to note from studies of Turing patterns are noted below. 

1. Specific kinetics are required for Turing patterns. Two component mechanisms, like the one de¬ 
scribed in Eq. can generate spatially heterogenous patterns. Whether or not systems are capable of 
generating spatial patterns through a Turing mechanism depends on the reaction kinetics. In partic¬ 
ular, we have outlined the conditions required in Eq. In general, these systems crucially have an 
activation-inhibition form (Eigure[^). Eurther details on such mechanisms can be found in ifT^ . 
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Figure 2: Region of parameter space in the Thomas system (Eq. that allows for Turing patterns for the 
ratio of diffusion constants D = Da/Du'- -D = 10 (blue), D = 100 (orange) and D = oo (green). The 
following parameters are fixed as follows: a = 1.5, p = 13, K = 0.125. As the D increases, the region 
over which Turing patterns can be observed increases. 

2. Diffusion can lead to instabilities under certain conditions. Given that the kinetics fulfill fhe con¬ 
ditions ouflined above, Turing insfabilifies can arise only if fhe rafio of diffusion consfanfs, D = 
Da/Du / 1. These spafial heferogeneifies arise especially when fhe diffusion consfanf of fhe ac- 
fivafor in fhe sysfem is smaller fhan fhaf or fhe inhibitor (called LALI, or ’’local acfivafion, laferal 
inhibifion” mechanisms). The concepf of diffusion-driven insfabilify was a groundbreaking concept 
when it was first proposed because diffusion was long-considered to be a stabilizing presence. The 
patterns that arise from these instabilities are called Turing patterns, and are thought to commonly 
arise in several biological systems. 

3.2 Cell polarization 

Polarization is fundamental to eukaryotic motility, morphogenesis, and division . Cell polarization is the 
process of reorganization of the cytoskeleton into distinct front and back regions in response to a stimulus 
. Recently, Mori et al. ifTTll presented the ‘wave pinning polarization’ or WPP model, a minimal model of 
cell polarization, that demonstrated both bistability and spatial patterning. In this system, the kinetics of the 
two species are given by: 

f{a, u) = u {ko + ^2 - 

g{a, u) = -f{a, u) = -u - koffO^ ■ 

After their initial proposal of the WPP model ifTTl . the authors considered several important extensions 
El [111131. However, a full characterization of the equilibrium behavior of the system as a function of the 
basal rate of activation /cq and the average amount of total protein p has not yet been adequately discussed. 
In this section, we outline such a characterization. The mathematical basis of such analyses can be found 
in the fundamental book by Strogatz 1251 . Computations were done in Mathematica and MATLAB; all 
associated scripts are available as a user-friendly GUI at www.ocf.berkeley.edu/~jnirody^^. 
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Figure 3: Equilibrium curves for each of the three steady states in the WPP model. For all of the following 
plots, we choose: 7 = 1, feoff = 1, K = 1 and vary the average amount of total protein p. The two stable 
steady states and are shown as solid blue and green lines, respectively; the unstable steady state 
is shown as a dashed line. For a range of feo, all three steady states exist and are real-valued; this region is 
shaded in red; this range increases with p, eventually resulting in an irreversible system response when it 
reaches feo = 0 . 


3.2.1 Well-mixed model 


As before, we begin by analyzing the reaction-only ODE model before considering the full reaction-diffusion 
system. The WPP model has three steady states, which we denote (a^^, n^^), (a^, n^), and ( 03 , To assess 
the stability of each of these steady states, we look at the Jacobian matrix of the system: 


J = 


2auK^') 
{a'^+K'^f ■ 
2auK^'y 


- feoff 
+ feoff 


ko + 1^ 

-ko-T^ 


(9) 


Example equilibrium curves for 7 = 1, feoff = 1, AT = 1 are shown for various values of average total 
protein p in Figure For the WPP model, and 03 are stable (shown as blue and green solid lines, 
respectively), while is unstable (shown as a black dashed line). For certain choices of parameters, both 
stable steady states exist; this is called the bistable regime. Within this region, either of the two steady states 
may be reached for the same set of kinetic parameters, depending on the initial conditions of the system 
(i.e., the ‘starting’ concentrations). 

Bistability is a common feature in biochemical reaction networks, particularly those containing positive 
feedback loops ll^ . Through bistability, positive feedback loops may allow for a sustained cellular response 
to a transient external stimulus lISTl . a central feature in cell polarization. To illustrate this, consider the basal 
rate parameter feo as a function of an external stimulus S (i.e., feo = k^S). Now, we can look to the plots 
in Figure as dose-response curves — the cell responds to external stimulus S by producing the activated 
protein A. 

As S (and consequently feo) is slowly increased, the concentration of A follows along the curve corre¬ 
sponding to 03 (blue) until it crosses the bistable region, after which the equilibrium value is the larger 03 
(green). If the stimulus is removed, and the level of S decreases, the higher-valued equilibrium is maintained 
within the bistable region; this behavior, where the dose-response relationship is in the form of a loop rather 
than a curve, is called hysteresis. If the bistable regime is large enough, in particular if it extends to fe = 0, 
an essentially irreversible response to transient stimuli may be elicited (as is seen for p = 0.31 in Figure]^. 

Figure shows the bistability region for the WPP model as the total average protein concentration p 
and the basal activation rate feo are varied. From the sloping shape of the bistable region, we can see that 
for sufficiently high values of total protein (for the set of parameters in Figure]^ p > 3), an ‘irreversible’ 
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Figure 4: The extent of the bistable region in {ko,p) space for the well-mixed model. In unshaded re¬ 
gions, only a single stable steady state (either ( 03 , or (a^,tt^)) exists. Within the shaded region, both 
steady states exist and are stable. The system may tend to either state, dependent on initial conditions. The 
parameter set used in the simulations by Mori et al. ifTTl is shown as a purple point. 


response may be generated: both steady states are stable even when the stimulus is removed, at ko = 0. The 
simulation parameters chosen by Mori et al. are shown as a purple dot. 

3.2.2 1-D spatial model 

Having characterized the bistable regime of the well-mixed model, we now turn our attention to the full 
reaction-diffusion system. The homogenous equilibria also serve as steady states for the spatially extended 
system and as before, we are interested in their stability throughout the parameter space. 

Recall that for the spatially homogenous system, we turned to the Jacobian to analyze the linear stability 
of these states under small perturbations from equilibrium. The eigenvalues of J are 

(/a -fu- {Da + Du)k^) ± \l{fa -fu- {Da + Du)k^)^ - 4 {DADuk^ + - Dufa)k^) 

( 10 ) 

As the expressions for the eigenvalues now contain the additional unknown k, we are now interested in a 
family of solutions, one for each wavenumber. The eigenvalues for the nonspatial system are achieved when 
fe = 0 , and are given by a~ = fa — fu and = 0 ; therefore any spatially homogenous perturbation will 
relax back to the spatially uniform steady state. We are instead interested in the emergence of heterogenous 
patterns that occur via perturbations within a finite range of critical wavenumbers 0 < kc < ^max- 

The region of ‘linear instability’ is highlighted in orange in Figurej^for A: = 0.2 pm~^. This corresponds 
to a pattern of wavelength w = ^ 30/rm. This length scale is intermediate among the motile eukaryotic 

cells that use Rho GTPases to generate polarity. In this region, one or both of the steady states, {a\,u\) and 
( 03 , U 3 ), lose stability with respect to a inhomogeneous perturbation as shown above (see Figure]^. 

We note that the region computed is for diffusion coefficients Djj = 10 and Da = 0.1 

pLrc?s~^. A similar parameter topology (with slightly larger regions of instability) was found in a reduced 
one-species model where infinite cytoplasmic diffusion was assumed . The disparity in diffusion coeffi¬ 
cients assumed in Figure [^presupposes the compartmentalization of the two species: a protein diffuses far 
more slowly on the membrane than in the cytosol (here, we assume the ratio of diffusion coefficients to be 
Ri 0.01 1201 ). In the following sections, we will assess the importance of this presumption. In addition to 
the region of linear instability, the parameter space for the full model can admit a surrounding region where 
front-like solutions are observed. In this region, a stalled wave can appear when the system is subjected to a 
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Figure 5: Parameter space topology for the full PDE model when Djj = 10/rm^s“^ and = 0.1/rm^s“^. 
The region of linear instability is shown shaded in orange for wavenumber A; = 0.2 This corresponds 

to a perturbation of length L = k, 30/um. Smaller values of k result in an expansion of the linear 
instability region; larger values of k result in the region shrinking. An additional domain is shown shaded 
in blue, in which front-like solutions are supported when given a sufficiently strong (or spatially graded) 
perturbation. The parameter choice made by Mori et al. (purple point) lies in this region. 


directed stimulus (e.g., a gradient) or if the domain exhibits some intrinsic polarity at f = 0. The parameters 
chosen by Mori et al. fall within this region; their simulations demonstrate that this ‘intrinsic polarity’ may 
arise via sufficiently noisy initial conditions iflTl . 

The boundaries of this regime are calculated by solving the Maxwell condition and finding fhe ranges 
of u fhaf admif a sfalled wave iflTl : 


Hb) = [ f{a,u)da = 0. (11) 

J a_ 

We fhen use fhe mass conservation condifion fo compute fhis region in (fco,p) space. While fhis is nol ana- 
lyfically feasible for fhe WPP model, numerically solving fhe infegral provides an accurafe characterization 
of fhis region. The mefhod ouflined here is nof specific fo fhe WPP model, or even fo models of cell polarity. 
To make such a phase-space analysis more accessible fo fhe biological communify, we provide an easy-fo 
use GUI fo allow readers fo analyze fhe sfabilify properfies of ofher sysfems of inferesf. 

3.2.3 Key take-away points 

The analysis above provides us significant insight into the behavior of the WPP model throughout the pa¬ 
rameter space (Figure]^. Using this, we can point out several notable properties of the model as they pertain 
to cell polarization. 

We note that the below properties are predicted using analyses performed on a ID spatial model. Ex¬ 
tension of the model to three, or even two, dimensions may (and likely do) result in different behaviors ||5l. 
However, an analytic treatment of higher-dimensional spatial models is rarely possible, and thus the com¬ 
prehensive analysis of a reduced model in one spatial dimension lays the foundation for simulation studies 
in higher dimensions. 

1. Importance of compartmentalization: The Rho GTPase family is large and varied, and is present 
in eukaryotes spanning from C. elegans to humans. However, one common feature of these proteins 
is compartmentalization-. the active form is bound to the membrane, while the inactive form diffuses 
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Figure 6: Examination of the rise of Turing patterns by loss of stability given critical wave numbers k. 

(A) A small segment of the parameter space is highlighted for the full spatial WPP model, showing regions 
where neither (region 1), one (regions 2 and 3), or both (region 4) of the equilibrium points become unstable. 

(B) An illustration of the loss of the linearly instability regime as D approaches 1. Plots show the magnitude 
of the real part of the rightmost eigenvalue for both equilibria within each of the regions highlighted in A. 
Plots are shown for A: G [—1,1]. When D = 10“^, corresponding to the localization of the active form to 
the membrane, a finite range of critical wavenumbers is observed; this range disappears when D = 1. 


in the cytoplasm ll2T]| . This feature has been shown to be important for cell polarization llOl fTTlfT^ . 
We illustrate the necessity of membrane localization by considering how the phase space of the WPP 
model in Figure [^changes if both species are contained in the cytoplasm. 

Qualitatively, the dependence of polarization on compartmentalization is relatively intuitive. As A 
and U constitute GTP- and GDP-bound versions of a single protein, their cytoplasmic diffusion rates 
are likely very similar. Given this, one would not expect the formation of any sort of regular pattern 
with no initial spatial structure. 

We show this quantitatively by defining the ratio of the diffusion rates D = Da/D u, and consid¬ 
ering the effect of this quantity on the regimes allowing spatially heterogenous solutions. For one 
membrane-bound and one cytoplasmic species, we take this ratio to be 0.01 Il 20 l . 

When D = 1, no spatially heterogenous patterns can be generated, and this region disappears alto¬ 
gether. However, as D decreases, spatial patterns are supported for a finite range of wavenumbers 
(Figure § iH. As the rate of diffusion of U increases, the maxima of a{k) move towards k = Q, 
eventually displacing the uniform steady state in the limit Du oo |[28l . This trend is not symmetric 
— as D increases from 1, there is no consequent extension of the multistable regime. This is similar 
to the formation of Turing patterns in local excitation, global inhibition models; for an overview of 
this brand of models with respect to cell polarization, we refer the reader to several excellent reviews 

Eiiiini. 

2. Benefit of being big: In addition to the importance of different diffusion rates between active and 
inactive forms, we can use the fact that the range of critical wave numbers is bounded above to 
consider the existence of a corresponding lower bound on the length of the cell L. 

The value of the maximum critical wavenumber /cmax for a feasible value of D is quite low (Figure]^. 
This suggests that smaller cells are less sensitive to polarization, while larger cells are able to respond 
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more robustly. Interestingly, this result has been observed experimentally: cells were found to become 
significantly more sensitized as they were flattened in a confined channel Q [TQ . 

3. Spontaneous polarization: Recall that fco can be written as a function of the concentration of some 
stimulus S: ko{S) = k^S. This formulation allows us to characterize parameter values which allow 
for spontaneous polarization in the absence of a directed external stimulus. 

Certain, but not all, cells are able to spontaneously self-polarize. In the WPP model, we can by 
considering the regimes crossed in Figure|^by the line ko = 0. We see that polarization in the absence 
of a stimulus is possible if the value of p is sufficiently high. This is consistent with experimental 
observations that the constitutive expression of Rho GTPases result in extension of randomly oriented 
lamellopodia and membrane ruffling 1 (221 . 

4. Polarization strategies: The parameter space topology for the WPP model contains two distinct 
regions that allow for non-homogenous equilibrium solutions. Because of the choice of parameters 
in Mori et al. mu, the system behavior in only one of these regions, corresponding to stalled-wave 
solutions (shaded in blue in Figure]^ was explored. The fact that these solutions were not initially 
found by a purely simulation-based study further emphasizes the need for comprehensive analytic 
treatment of biological models. 

In general, Turing patterns form more easily (i.e., in response to far smaller perturbations) than pat¬ 
terns formed by a wave-pinning mechanism. However, they occur on a far slower timescale ElEl. 
The existence of a Turing-like instability regime in addition to a region which admits stalled-wave 
solutions presents cells with multiple strategies for polarization. 


4 Conclusions 

Mathematical modeling has served to complement experiment to further our understanding of biological 
systems. Recent advances in knowledge of these systems have allowed us to construct more and more 
complex models, but our computational abilities have not caught up. In particular, our knowledge of the cell 
has made us increasingly aware of the structured nature of the cellular environment. 

Despite this knowledge, biological models often neglect this structure and consider biological networks 
in a homogenous spatial environment. These purely kinetic models are computationally more tractable, but 
can fail to capture the biological reality. However, incorporating spatial structure into these mathematical 
models can make them difficult to analyze using current tools. 

For systems involving spatial patterning, such models are often of the reaction-diffusion type. As the 
behavior of these systems can vary widely with parameter choice, simulation studies are often unsatisfying 
in that they only provide a small peek into the full capabilities of these models. Above, we have outlined the 
steps for a comprehensive analytic treatment of spatial structure arising from reaction-diffusion models. 

Additionally, we described examples from two distinct patterning mechanisms: diffusion-driven insta¬ 
bility (also known as the Turing mechanism) and wave-pinning. We have highlighted the usefulness of such 
a protocol, in that it illuminates several features of models which are difficult to observe via simulation 
alone. 
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